Ocean warming and acidification adjust inter- and intra-specific variability in the functional trait expression of polar invertebrates

Climate change is known to affect the distribution and composition of species, but concomitant alterations to functionally important aspects of behaviour and species-environment relations are poorly constrained. Here, we examine the ecosystem ramifications of changes in sediment-dwelling invertebrate bioturbation behaviour—a key process mediating nutrient cycling—associated with near-future environmental conditions (+ 1.5 °C, 550 ppm [pCO2]) for species from polar regions experiencing rapid rates of climate change. We find that responses to warming and acidification vary between species and lead to a reduction in intra-specific variability in behavioural trait expression that adjusts the magnitude and direction of nutrient concentrations. Our analyses also indicate that species behaviour is not predetermined, but can be dependent on local variations in environmental history that set population capacities for phenotypic plasticity. We provide evidence that certain, but subtle, aspects of inter- and intra-specific variation in behavioural trait expression, rather than the presence or proportional representation of species per se, is an important and under-appreciated determinant of benthic biogeochemical responses to climate change. Such changes in species behaviour may act as an early warning for impending ecological transitions associated with progressive climate forcing.

maintain 38 , reduce 35 or enhance [39][40][41] functioning, making it difficult to generalise species contributions to alterations in ecosystem properties.Disentangling these effects is frustrated by the fact that changes in behaviour are also accompanied by compensatory responses 42,43 that affect dominance patterns 44,45 , and other factors, which can partially, or wholly, offset functional responses to forcing 46 .Nevertheless, field observations show that a shift in the type and amount of faunal activity can lead to environmental transitions 3 that exert a disproportionate influence on ecosystem properties and functioning over and above the effects caused by changes in species diversity 47,48 and composition 45,49 .It is important to note, however, that although flexible behavioural strategies can improve short-term fitness 50,51 , any associated functional consequences 52,53 may not materialise until much later and can be hard to distinguish from other temporal changes in the system 54 .
We anticipated that changes in species behaviour will be more pronounced in regions of fast paced climate change 3,55 , as genetic and other coping mechanisms are less likely to be enacted in time.We speculated, given the closure of dispersal and adaptation as viable options, that adjustments to individual behaviour would dominate species responses to change 56 at higher latitudes.Here, using sediment-dwelling invertebrate species obtained from the Arctic Barents Sea (the bivalve Astarte crenata, sea star Ctenodiscus crispatus and polychaete Cistenides hyperborea) and Antarctic Peninsula (the protobranch Aequiyoldia eightsi and bivalve Laternula elliptica), two areas currently experiencing amplified climate change 57,58 , we explore the combined effects of near-term ocean warming (+ 1.5 °C) and elevated levels of atmospheric carbon dioxide (550 ppm [CO 2 ]) on aspects of species behaviour known to influence biogeochemical cycling.As we anticipate that the direction and magnitude of change in behaviour will diverge between species 4,59,60 , we also include individuals of Astarte crenata and Ctenodiscus crispatus from two locations within the Barents Sea that contrast in temperature and sea ice dynamics; here, our expectation is that individual species responses will be in line with previous observations 3 , but will be more pronounced when species are from locations experiencing narrower environmental variation.We use these data to demonstrate the importance of behavioural change and compensatory mechanisms, including numeric and/ or biomass increases and performance enhancement 42,43 , in moderating how benthic environments respond to external forcing.We show, for five species of polar benthic invertebrates, that the ability to modify behaviour in the face of climatic forcing does not guarantee that species contributions will remain unchanged.Our findings emphasise the importance of context-dependency and have implications for the functional contributions of populations facing climate change, their capacity to adapt in the face of further environmental transitions, and suggest that the onset of phenotypic expression may serve as an early warning for impending ecological change.

Results
We find evidence that individual movement and burial behaviour, sediment particle reworking activity, burrow ventilation activity, and associated nutrient concentrations at the sediment-water interface, can be dependent on environmental condition (ambient climate treatment vs future climate treatment of + 1.5 °C and 550 ppm [CO 2 ]), location, and species identity (Supplementary Models S1 to S29).However, observed effects seldom form full factorial interactions between the three dependent variables (8 of 29 models).Despite observing mortalities in the bivalve Astarte crenata (2 individuals, 1 from each climate), the sea star Ctenodiscus crispatus (4 individuals, 3 ambient and 1 future climate), and the polychaete Cistenides hyperborea (1 ambient climate), it was possible to relate our response variables in ecosystem process (sediment particle reworking: surface boundary roughness, median mixed depth and maximum mixed depth; burrow ventilation activity) and functioning (nutrient concentrations: ammonium, nitrite, nitrate and phosphate) to species behaviour (individual movement: response time; burial behaviour: burial time) in all aquaria.We find no evidence that differences in mortality (assessed using total biomass as a random effect) affects our response variables.

Discussion
Our findings demonstrate that conditions representative of anticipated near-future climate change can lead to fundamental shifts in functionally important aspects of sediment-dwelling invertebrate behaviour.These effects can be substantive; here we observed a doubling of burial rate, deepening of particle mixing and a change in the magnitude and direction of biogeochemical dynamics that are sufficient to change the functional role of a species (A.eightsi 36 ).This observation is important, because alterations in individual functional capacity that are distinct from functional shifts caused by changes in community composition and/or novel environmental conditions are common 3,61 , and likely result from changes in the strength and nature of a portfolio of sublethal responses, including species interactions 62,63 , compensatory mechanisms 41,42 and/or other subtle phenotypic responses 54,64 .Changes in macronutrient cycling under climate forcing is not trivial to detect 65 , however, and may be masked by the pH buffering effects of [CO 2 ] driven alkalinity changes 66 on microbial mediated pathways of nutrient recycling.
The behavioural changes with warming and acidification observed here may be even more important ecologically in polar regions than they would be at lower latitudes.Seasonality results in many species entering low energy and activity states similar to aestivation in winter that can last several months 67,68 ; in this study, as in L. elliptica 69 , although juvenile A. eightsi growth is known to be similar across summer and winter 70 .Therefore, in the presence of species that respond to seasonal cues, greater levels of species activity, leading to greater microbial and nutrient remobilisation from sediments 32,33 , may occur for longer in polar regions as the summer season extends under climate change 71 .If widespread, it follows that there will be positive ramifications for phytoplankton productivity over the long term 1,3 .Although this is not the only mechanism underpinning nutrient provision for productivity, we speculate that outcomes associated with benthic responses to climate change could include changes in the phenology of the initiation of productivity and early intensity of phytoplankton growth 72 , with downstream impacts for primary and secondary consumers.
Whilst the effects of a near-future climate in our experiments were comparatively weaker than the effects of species identity and location, consistent with theoretical expectations 73,74 , we did note a reduction in intraspecific variation that reflected changes in environmental context and location 37 .This can be very important for maintaining populations 75 , enabling adaptation to changing environmental conditions 76 and stability in ecosystem functioning 77 .However, whilst sublethal responses may enable species to persist in, or for longer, under novel circumstances, other phenotypic costs may constrain or inhibit an individual's ability to adjust further 78,79 .Hence, reductions in intra-specific variation may serve as an early warning for impending ecological transitions associated with progressive forcing and potentially inform more timely management actions, reinforcing the need for continual monitoring of faunal activity and the ecological constraints that modify functionally important aspects of species behaviour 80 .
The variation in intra-specific behaviour observed here under enhanced warming and [CO 2 ] is consistent with other behavioural studies 81 , physiological responses in polar benthic species 21 and incorporates regional contextualisation 13 .Whilst our study was not explicitly designed to examine species range shifts or gradients of environmental change, an important feature of our sampling design was that our locations were positioned to the north and south of the oceanographic polar front, contrasting in benthic biogeography 82 , bioturbation activity and functioning 3 .Hence, we were able to show that individuals predisposed to a wider inter-annual thermal range exhibit a more reserved behavioural response to change than those inhabiting a narrower thermal range 83 .Thus, plasticity in response mirrors the level of local environmental fluctuation 84 .Whilst spatial associations between environmental temperature range and physiological thermal tolerances are not atypical in ectothermic species 13,85,86 , this does mean that high latitude populations may be at greater risk of local extinction over the long term.As thermal tolerance narrows with decreasing seasonality in temperature towards the poles 16,87 , and will likely be further constrained with ocean warming 88 , populations already at or approaching the edge of their thermal limits will most likely have less scope to compensate and adapt to change 89 .Indeed, changes in species composition and abundance are well documented across areas of environmental transition 3 and show similar patterns of functional change, as observed here.Temperature-driven responses are, however, typically complicated by interactions with other abiotic drivers 74 and are likely to lead to both amplified and dampened effects in spatially stochastic ecosystems 90 .Yet, previous investigations have predominantly focused on spatial distributions of species turnover 64 , functional diversity 91,92 and redundancy 93 , rather than characterising intraspecific variability of species-environment interactions.The latter can be a more important driver of the short-term functional response of communities than changes in species composition, dominance, or richness 94,95 .For example, the shallower burrowing activity of invertebrates held under more acidified conditions 96 allows species to evade the physiological effects of decreasing pH, but simultaneous burrowing and ventilatory 40 responses to warming to maintain environmental continuity may negate the need for such avoidance behaviour 97 .We observed similar changes across multiple aspects of functionally important behaviour that may have led to non-additive effects on net functioning that were not possible to distinguish.Nevertheless, the cumulative effect of such short-term behavioural responses is likely to be decisive for the composition 28 , population dynamics 98 , connectivity 99 and functioning 100 of benthic communities that will be moderated by seasonal timing 54 and local circumstance 13,36 , including interannual variability 3 .
Quantitative information on the functional role of individual species is rare for both polar regions 101 , yet understanding, and accounting for, species responses to climate change is fundamental to improving the likelihood of determining the most realistic ecosystem future 102 .We contend that this task will be frustrated by context-dependent variation in both intra-and inter-specific responses to forcing that are not readily captured using fixed trait modalities 35,103 .Where the overall outcome of species responses remains largely unresolved, reductions in the variation of conspecific responses 95,104 may form a viable alternative for some predictive models.Furthermore, our findings lend support to the view that location-dependent variation in behavioural www.nature.com/scientificreports/responses can be attributed to localised thermal plasticity driven by exposure to divergent temperature seasonality trends 8,84,105 .Inter-and intra-specific variations in vulnerability, effect-and-response traits 79 and interactions between species 106,107 can facilitate functional redundancy and/or post-change compensations 42,43 .A mechanistic approach that explicitly tests suspected abiotic and biotic signals is necessary for establishing patterns of response 108 across multiple levels of biological organisation 109,110 , enabling the generation of more robust projections of the most likely functional consequences of change.

Fauna and sediment collection
We obtained individuals of the bivalve Astarte crenata, sea star Ctenodiscus crispatus and polychaete Cistenides hyperborea from replicate SMBA (Scottish Marine Biological Association, 50 × 50 cm) box cores, and 15 min Agassiz trawls in the Barents Sea (stations B13, 74.3° N, 30.0°E; B16, 80.3° N, 30.0°E, 263-375 m depth; JCR18006, RSS James Clark Ross, Supplementary Fig. S1a, Table S1) in July 2019.Individuals of the protobranch Aequiyoldia eightsi and bivalve Laternula elliptica were collected by SCUBA divers at Rothera Point, Adelaide Island, West Antarctic Peninsula (67.3° S, 68.1° W, 10-20 m depth, Supplementary Fig. S1b) in March-April 2019.We obtained surficial sediment (< 5 cm depth) from SMBA box cores at the Barents Sea stations B13, B14 and B16 (Supplementary Table S1) for the Arctic species, and from the intertidal mud flats of the Hamble, UK (50.9°N, 1.3° W) for the Antarctic species.Each sediment was sieved (500 µm) within a seawater bath to retain the fine fraction and to remove macrofauna and debris.Sediment particle size (Supplementary Fig. S2) was determined using a Malvern Mastersizer 2000 He-Ne LASER diffraction sizer.Mean particle size, sorting, skewness and kurtosis were quantified using GRADISTAT 111 .Loss on ignition was used to determine sediment organic matter content (%).

Experimental design and set-up
Sediment (mean ± s.e., n = 38: particle size = 60.30± 3.91 µm, organic matter content = 5.502 ± 0.212%; Supplementary Table S2) and species were distributed across 42 clear acrylic aquaria (internal LWH: 12 × 12 × 33 cm, 3 replicates treatment −1 : species × location × climate scenario; Supplementary Table S1), designed to accommodate representative field densities (Arctic species, 2 ind.aquarium −1 ; Antarctic species, 1 ind.Aquarium −1 ; ( 112 ; Supplementary Table S4) and the size and burrowing requirements of each species (sediment depth: A. crenata, C. crispatus & C. hyperborea, 16 cm; A. eightsi, 12 cm; L. elliptica, 19 cm 113,114 ).Aquaria were randomly placed within one of two insulated seawater reservoirs ( 3 , Supplementary Fig. S3) in the Biodiversity and Ecosystem Futures Facility, University of Southampton (UK).All aquaria were filled with seawater (salinity 33, 10 µm sand filtered, UV sterilized) to ~ 12 cm above the sediment-water interface and maintained in the dark.After a transitionary period to aquarium conditions (21  ).This period of time was sufficient for the establishment of microniche formation 116 and vertical biogeochemical gradients indicated by colour change 117 to form in the sediment.Partial seawater exchanges (weekly, 50% volume) prevented accumulation of excess food and nutrients.Measurements in behaviour, ecosystem process and functioning were taken at the end of the experimental period.

Seawater carbonate chemistry, temperature, and salinity
Atmospheric [CO 2 ] (Supplementary Fig. S4) was controlled using a custom-made CO 2 -air mixing system which continually maintained and monitored [CO 2 ] in the air mixture supplied to each individual experimental core using infrared analysers (LI-COR LI-840A) 54 .This approach facilitates natural variability within the carbonate system 118

Behavioural response of individuals
Behaviour of C. crispatus, C. hyperborea and A. eightsi were quantified using measurements of movement and burial behaviour at the sediment surface.Individuals (morphology, ± 0.01 mm; blotted wet weight, ± 0.001 g, Supplementary Table S5) were placed in separate treatment-acclimatised viewing trays containing sediment (depth 5 cm) overlain with sea water (depth 3 cm) and viewed (≤ 60 min) with a benchtop video camera (Logitech C920 HD Pro, 1080p; Supplementary Fig. S7).The time taken to initiate movement (response time, s) and to complete burial (burial time, s) was recorded (3 frame s −1 , SkyStudioPro) and analysed frame-by-frame (VLC Media Player).We incorporated biomass as a random factor in the statistical analysis to account for any intra-specific variation in size.
Vol:.( 1234567890 -S12), and the distribution of luminophores was analysed using ImageJ (version 1.46r 120 ).From these profile data (Supplementary Fig. S13), we calculated the mean ( f-SPI L mean , time dependent indication of mixing), median ( f-SPI L med , typical short-term depth of mixing) and maximum ( f-SPI L max , maximum extent of mixing) mixed depth of particle redistribution.Given the shape of the vertical distribution of luminophores (non-continuous), f-SPI L mean was an unsuitable descriptor of the distribution profile and not considered for statistical analysis.The rugosity of the sediment-water interface (upper-lower limit = surface boundary roughness, SBR) provides an indication of surficial activity.Ventilation behaviour 101 of all five species was estimated from absolute changes in the concentration of sodium bromide [NaBr] 54 .Dissolved [NaBr] was standardised across all aquaria (mean starting concentration = 1353.816± 317.264 mg L −1 ) and [NaBr] was determined using a Tecator flow injection auto-analyser (FIA Star 5010 series).Negative values of [NaBr] (∆[Br − ] mg L −1 ) over an 8-h period indicate increased infaunal ventilatory activity.
As faunal activity mediates nutrient concentrations, we determined water column [NH 4 -N], [NO 3 -N], [NO 2 -N] and [PO 4 -P] (µmol L −1 , ~ 10 ml, filtered 0.45 μm NALGENE nylon matrix) for all five species once a month (Supplementary Fig. S14) using a QuAAtro 39 auto-analyser (SEAL Analytical) as a measure of ecosystem functioning.As nutrient concentrations will also reflect differences in the volume of sediment between species treatments, we calculated the log response ratio (lnRR = ln[conc before /conc after ] 121 ), an effect size that quantifies proportionate change.As patterns of [NO x -N] are reciprocal to those of [NH 4 -N] but indicate beneficial biogeochemical processes (e.g.denitrification), lnRR values for [NO 2 -N] and [NO 3 -N] were multiplied by −1 to align the direction of ecosystem functioning.

Statistical analysis
Analysis of Variance (ANOVA) models were developed for each dependent variable (movement and burial behaviour: response time, burial time; ecosystem process: SBR, f-SPI L median, f-SPI L max , ∆[Br − ]; ecosystem functioning: [NH 4 -N], [NO 3 -N], [NO 2 -N], [PO 4 -P]).For A. crenata and C. crispatus, we determined the effects of the independent variables; environmental condition (2 levels: ambient, future), location (2 levels: stations B13 and B16; Supplementary Fig. S1a), species identity (2 levels), and their interactions, whilst for A. eightsi and L. elliptica, we determined the effects, alone and in combination, of the independent variables environmental condition (2 levels) and species identity (2 levels).As C. hyperborea was found at a single station, we determined only the effects of environmental condition (2 levels).Intra-specific variability within treatment levels was determined using the coefficient of variation.
Model assumptions were visually assessed using standardised residuals vs fitted values plots, Q-Q plots, and Cook's distance 122 .Where there was a violation of homogeneity of variance, we used a varIdent variance-covariance structure and generalised least-squares (GLS) estimation 123,124 to allow residual spread to vary amongst groups.We determined the optimal fixed-effects structure using backward selection informed by Akaike Information Criteria (AIC) and inspection of model residual patterns.For the GLS analysis, we determined the optimal variance-covariance structure using restricted maximum-likelihood (REML) estimation by comparing the initial ANOVA model without variance structure to equivalent GLS models incorporating specific variance terms.These models were compared against the initial ANOVA model using AIC informed by visualisation of model residuals.We determined the optimal fixed structure of the most suitable model by applying backward selection using the likelihood ratio test with maximum-likelihood (ML) estimation 122,124 .For ANOVA models with interactions, we calculated the effect size (ω 2125 ) of each independent variable in R 126 using the effectsize package 127 .For GLS models with interactions, we determined the relative importance of each independent variable by comparing the minimal adequate model with a model with the independent variable of interest, and all its interactions, removed using likelihood ratio (L-ratio) in the nlme package 123 .Details of initial and minimal adequate models (Model S1 to S29) are provided in electronic supplementary material.